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Abstract 

Prognostics of large composite structures is a topic of increasing interest in the field of structural health monitoring for aerospace, 
civil, and mechanical systems. Along with recent advancements in real-time structural health data acquisition and processing for 
damage detection and characterization, model-based stochastic methods for life prediction are showing promising results in the 
literature. Among various model-based approaches, particle-filtering algorithms are particularly capable in coping with 
uncertainties associated with the process. These include uncertainties about information on the damage extent and the inherent 
uncertainties of the damage propagation process. Some efforts have shown successful applications of particle filtering-based 
frameworks for predicting the matrix crack evolution and structural stiffness degradation caused by repetitive fatigue loads. 
Effects of other damage modes such as delamination, however, are not incorporated in these works. It is well established that 
delamination and matrix cracks not only co-exist in most laminate structures during the fatigue degradation process but also 
affect each other’s progression. Furthermore, delamination significantly alters the stress-state in the laminates and accelerates the 
material degradation leading to catastrophic failure. Therefore, the work presented herein proposes a particle filtering-based 
framework for predicting a structure’s remaining useful life with consideration of multiple co-existing damage-mechanisms. The 
framework uses an energy-based model from the composite modeling literature. The multiple damage-mode model has been 
shown to suitably estimate the energy release rate of cross-ply laminates as affected by matrix cracks and delamination modes. 
The model is also able to estimate the reduction in stiffness of the damaged laminate. This information is then used in the 
algorithms for life prediction capabilities. First, a brief summary of the energy-based damage model is provided. Then, the paper 
describes how the model is embedded within the prognostic framework and how the prognostics performance is assessed using 
observations from run-to-failure experiments. 
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1. Introduction 


Recent advancements in real-time structural health monitoring (SHM) methods enable an any-time in-situ 
assessment of damaged and aging structures’ condition and allow collecting data on the progressive, inevitable 
damage growth due to operating and contingent loads. SHM technologies are increasingly used for fatigue-induced 
degradation monitoring of composite structures, where the damage may be hidden in the internal layers and barely 
visible on the outer surfaces. The availability of information on damage type and damage extent would facilitate also 
the prediction of the remaining useful life (RUL) of the structure through the estimation of the fatigue damage 
progression. Addressing the RUL prediction by real-time methods can revolutionize the maintenance policies of 
aeronautical, mechanical and civil industry, moving from the current damage tolerance approach to condition-based 
or predictive maintenance strategies. 


However, real-time prediction of RUL of composite materials is a challenging task that must factor in the 
coexistence of multiple damage mechanisms or multiple damage-modes (MDMs). These damage modes interact 
with one another and might also generate new damage modes. A typical example is the interaction caused by matrix 
cracks, delamination and buckling. Matrix cracks can induce local delamination, which can become global 
delamination, and the delaminated layers can fail because of buckling in case of negative load ratios (i.e., when the 
load falls below zero) [1]. In addition to the coexistence of MDMs, the RUL prediction is inherently affected by 
several sources of uncertainty. Fatigue of materials is uncertain in nature, since it is driven by inclusions and 
impurities caused by the manufacturing process and complex physical nano- and micro-scale phenomena not 
accounted for in common engineering models. Where damage is measured with using automated SHM tools, there 
are additional uncertainties about current damage location and extent that further complicate the prediction of the 
damage growth. This makes stochastic approaches a logical choice for real-time RUL prediction. 


While highly desirable, real-time damage prognosis of composite laminates has so far been only sparsely explored. 
It requires a methodology to predict the damage evolution that merges stochastic approaches and real-time 
diagnostic information in the prognostic stageto update the RUL prediction. Recently, Bayesian filters have shown 
promising results in predicting the evolution of matrix cracks and consequent stiffness degradation of cross-ply 
laminates by including SHM data in the prognostic stage [2]. However, the effect of co-existing damage 
mechanisms was not incorporated in the damage progression model, and the resulting stiffness degradation was only 
linked to the presence of matrix cracks. 


The present work reported herein follows the methodology proposed in [2] and extends the Bayesian framework by 
including an energy-based MDM model. This model was recently investigated in [3], where the authors showed 
successful predictions of the stiffness degradation caused by matrix cracks and delamination in cross-ply fiber- 
reinforced laminates. In this study here, the MDM model is embedded in a Bayesian filtering algorithm commonly 
referred to as particle filtering [4], which aims at simultaneously monitoring three degradation processes: matrix 
crack density, delamination and stiffness reduction. The energy-based MDM model enables the estimation of the 
interacting damage growth rates, and the embedding of the mechanical model into the Bayesian framework allows 
the probabilistic estimation of the RUL conditioned upon the available diagnostic data. The developed model-based 
prognostic framework is assessed against tension-tension fatigue damage progression data on a carbon fiber- 
reinforced polymer (CFRP) laminate. The growing damages were observed through X-ray radiographies and their 
extents are provided to the algorithm sequentially, thus simulating a real-time condition where a SHM system 
provides regular information on the damage extent as time passes by. 


The rest of the paper organizes as follows: Section 2 introduces the MDM model and the equations to predict the 
damage progression, Section 3 shows the probabilistic framework based on particle filtering and the basic 
probability density functions (PDFs) to implement the algorithm and Section 4 shows the application on CFRP 
damage growth data. Section 5 draws some conclusions. 


2. Modeling of concurrent damage progression 


Most of the models for damage growth prediction in composites resort on finite element methods because of the 
complexity of the damage mechanisms and their interactions. Though, the computational costs of finite elements 
commonly prevent the applicability in stochastic algorithms for real-time applications. Analytical damage 
progression models are therefore preferred in real-time prognostic scenarios: the model complexity can be scaled 
down to perform a fast estimation of the damage progression using simplified formulations of the stress state [3]. In 
addition, the tuning of analytical models is usually easier because of the limited number of model parameters. 


The approach proposed here, which follows [2], resorts on a strain energy release rate (SERR) model, since the 
stress intensity factor, widely adopted in metal fatigue, loses its usefulness when several cracks and different 
damage mechanisms affect the material. The formulation of the SERR includes the effect of matrix cracks and 
delamination, G = G(p,D), which are the two most common damage mechanisms affecting fiber-reinforced 
laminates. The model provides also an estimation of the Young’s modulus of the damaged laminate E,,, enabling the 


monitoring of the stiffness degradation. 


The SERR range AG, caused by fatigue loads, is used as input to power laws for the estimation of the damage 
growth rates, intended as damage growth per load cycle [5]. These power laws are commonly named modified 
Paris’ laws, given the similarity with the Paris’ law used for metallic alloys. The growth rates of matrix crack 
density p and delamination D are therefore expressed using equations (1) and (2), respectively. 
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y= A(AG(p, D)) (1) 
dD B 
Ww B(AG(p,D)) (2) 


The closed form solution of (1) and (2) rarely exists, since G is mostly a highly nonlinear function of the amount of 
damage. So, the alternative is the estimation of the damage progression using a linear damage accumulation rule (3)- 
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Where AN is the number of load cycles between two discrete time-steps, k-/ and k. Since G of composite laminates 
is particularly high during the first stage of the fatigue life, large time steps (i.e., large AN) can generate inaccurate 
solutions. Therefore, AN is usually kept equal to 1 to simulate the damage growth correctly. Once the amount of 
damage has been calculated, the ratio of the current elastic modulus E,, (depending on p and D) and the initial elastic 
modulus E;, determines the normalized remaining stiffness of the laminate, (5). 
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The work in [3] investigated the capability of both matrix crack-induced delamination models and edge 
delamination models in describing the SERR as a function of matrix cracks and delamination and in estimating the 
laminate remaining stiffness as well. The study pointed out as the edge delamination model proposed by Zhang, 
Soutis and Fan [6] appeared as the best in describing the stiffness reduction of a notched cross-ply laminate already 
utilized in [2]. Zhang et al. analyzed the pioneering work of O’Brien [7], introducing the effective elastic modulus of 
a partially-delaminated laminate, and they enhanced this model including the effect of matrix cracks in the laminated 
region. 


The model requires some simplifying assumptions on the type and shape of delamination, crack pattern and crack 
location. Specifically, this model is valid for symmetric balanced laminates under tensile loading and matrix cracks 
spanning all the width of the 90° plies. Figure 1 shows a graphical representation of the model, using a [0,/90n]; 
stacking sequence. Delamination D supposes to span all the height of the laminate (x-direction) and grow along the 
width, transversally to the applied load (y-direction). It is expressed in meters, m. The matrix crack density p is 
measured as number of cracks per unit length, #/m. 
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Figure 1. Zhang's model for a partially-delaminated cross-ply laminate. 


In this case, the outer sub-laminates are aligned with the direction of the remote tensile stress 09 indicated by the 
thick arrows (0° sub-laminates). The inner sub-laminate is perpendicular to the load direction (90° sub-laminate) and 
is affected by fatigue-induced matrix cracks. The laminate is split in three parts, two of them are equal to one 
another. So, the symmetric geometry generates two regions with different features. Region I is the laminated region 
where delamination is absent and the matrix cracks in the 90° plies reduce the sub-laminate stiffness. Region II is 
the fully-delaminated region where the 0° sub-laminates disconnect from the 90° sub-laminate. The two regions 
have different longitudinal elastic moduli named E,.; and E,,;;, respectively. The effective elastic modulus of the 
laminate is the weighted average of the two elastic moduli (6). 
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Ey, 1, depends on the number of delaminated interfaces. If both the interfaces between the 0° and 90° sub-laminates 
are delaminated (as has been shown in Figure 1), E,;; can be easily calculated from the laminate theory [1]. The 
matrix crack density p influences F,;, and this influence has been modeled in [6] by means of in-situ damage 
effective functions (IDEFs), A = A(p), which reduce the stiffness matrix Q of the 90° plies, (7). 


Q110 Q12,0 0 Qi2,0/Q22,0 Q12,0 0 A22(p) 
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The IDEFs Az2(p) and Agg(p) were determined by Zhang and coauthors in previous works [8]. They are based on a 
modified two-dimensional shear lag analysis and they describe the in-situ constraint conditions of the 90° plies. The 
equations to derive the IDEFs are not reported here for the sake of brevity. The interested reader is referred to the 
original papers [6]-[8] for finer details. Eventually, the SERR is calculated using (8), which already neglects any 
thermal effect. 


G(p,D) = (e(6,D))° * (Exi(9) ~ Exun) (8) 


Where e(p, D) = 09/E, (p, D). As explained above, the SERR G is embedded into the damage growth rate and the 
longitudinal elastic modulus of the damaged laminate E,(p,D) is used to estimate the remaining stiffness. It is 
worth noting that the SERR range AG has been defined using (9) and (10). In this way, the mean load effect, which 
can alter the similitude criterion, is neglected [9]. 


AG = (a Gnax ~ Gin )’ (9) 
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The SERR range in (9) feeds equation (1) and (2) to calculate the damage growth rates, and the model underneath 
particle filtering is then composed by equations (3)-(5). 


(10) 


3. Particle filtering-based prognosis 


This section describes the prognostic framework that composes of the MDM model in Section 2 and particle 
filtering. Specifically, a bootstrap sequential importance resampling algorithm with systematic resampling has been 
chosen to build the Bayesian framework [10]. Furthermore, a sub-algorithm has been used to update the model 
parameter during run-time. The fundamental equations of the algorithm are summarized and directly tailored for the 
fatigue damage prognosis problem below. 


Let us consider the system’s state x governed by a dynamic state-space (DSS) model (11). 


Xe = f (Xp~1, 9, Up_1, Ox-1) 
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The system’s state vector contains the three damage modes x, = [Px,Dx,5;,]7, while the measurement vector 


contains the observations of the true (unknown) system’s state z, = [Pr Deo Sk] - The evolution equation f 
describes the system’s state dynamics and is driven by equations (3)-(5). Instead, the observation equation g links 
the measures with the damage state. The uncertainties affecting the system’s state dynamics and the measurement 
system are embedded in the DSS using random processes called process noise w,, and measurement noise 7);, 
respectively. Since both the system’s state vector and the measurement vector are three-dimensional, the noises have 
been split in three independent contributions. Following the discussion in [11], the errors affecting the matrix crack 
density and delamination models are defined as log-Normal random processes e®, while the error of the stiffness 
degradation model is an unbiased Gaussian random process (12). The measurement noises are modeled as 
independent, unbiased Gaussian processes (13). 
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Equation (12) shows the dependence of the damage growth rates on the model parameter vector @ and the input 
vector u. The latter is the far-field stress range that drives AG, so u—> u = Ado = Omax — Omin. AS already 
presented in [11], the use of a log-Normal random process with specific relation between mean and variance 
produces an unbiased evolution equation. The process noise w follows equation (14), while the measurement noise 
is unbiased Gaussian 9 ~ MVN(0, 2): 
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The process and measurement noises have been modeled as stationary random processes. Therefore, their 
dependence from the time step k has been neglected. The model parameter vector @ can be updated during run time 
to improve the prediction performance of the algorithm. Here, the artificial dynamics sub-algorithm has been used 
for its simplicity and effectiveness [12]. It introduces a perturbation of the model parameter values using a random 
disturbance y, usually defined as an unbiased Gaussian process (15). 


Oy = Ox-a + VR-1 (15) 


Where y, ~ NV (0, Bae) The covariance matrix Zy;, must decrease as time passes by to guarantee the convergence 
of the algorithm [12]. In this application, the variance follows equation (16), chosen empirically according to the 
authors’ experience. 


1 
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Where %,, is the initial covariance matrix and k is the time step. The model parameter vector contains the empirical 
parameters of the modified Paris’ laws in equations (1) and (2), 8 = [log A, a, log B, B]’. The two parameters A and 
B have been embedded using their logarithmic form, since they are log-Normally distributed [13]. 


The sequential importance resampling method allows the approximation of the conditional PDF of the system’s state 
and model parameter given the observations, p(x, 0, |Zo.x), using N, weighted samples (17). 
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Where the pairs fae 0)" are Monte Carlo samples of the system’s state and model parameters, which are 
i=1 
) 


weighted by the weights wit di = 1,...,N,. Following the bootstrap particle filtering theory, the samples rs and 
oes are generated through the evolution equation and the artificial dynamics equation, (18). 
oe (i) 
x ,Uz_1,@ 
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Where yo? is a sample from NV (0, Brie) and w© is a sample from MVN(u,,,,,)- It should be noted that the 
subscript k on the model parameter vector refers to the pu that the samples are from the k-th posterior estimation, 
because the true @ is not time-varying. The weights wy? depend on the likelihood of the observation given the 


sample p(zx |x? ys and are then normalized to sum up to 1, (19). 
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The prognostic stage is carried out by propagating the samples of the system’s state in the future using the evolution 
equation f(-). At each time step, the system’s state sample x is altered by a sample of the model error w©, which 
simulates the unpredictable fluctuations of the damage growth rates caused by small-scale phenomena neglected in 
the model. The prognosis is analytically expressed through the p-steps ahead prediction equation (20), [14]. 
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Where p(x;|xj-1) is the transition density function, which comes from the probabilistic form of the evolution 
equation f (-). Once the posterior distribution of the system’s state has been computed, the samples ao are projected 
in the future using the transition density function (i.e., the evolution equation). The concept behind the prediction 
equation (20) can be further extended to calculate the number of fatigue load cycles to reach a pre-determined 
critical threshold, xcr. The samples of the system’s state are propagated using the transition density function until 
they reach the threshold, i.e., 8 = Xcr. The number of fatigue load cycles corresponding to the time step k + l is 
defined as the end of life of the sample, Nee Thus, the RUL of the i-th sample is the difference between Ni? and 


the number of load cycles of the current time step, RUL® = NO — N,. Once all the samples have reached the 


; Ns 
critical state, the sample pairs {RULP, we? enable the estimation of the RUL distribution, (21). 
i=1 


p(RUL,|Zo:¢) ~ Es, we?.5(RULW — RUL,) 
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Here, the critical threshold has been defined as a limit damage state xcr = [Pcr, Der, Scrl’ that cannot be exceed to 
guarantee the safety of the structure. Once the i-th sample x0 reaches one of the limits (either pcr, Der or Scr), the 


sample’s propagation stops and the related RUL® is recorded to compute the RUL distribution. 


4. Application to real fatigue-induced damages in CFRP laminates 


The algorithm in Section 3 is applied to fatigue damage progression data publicly available in the NASA prognostics 
data repository [15]. The picked experimental data refers to the notched cross-ply coupon LJS// with dog-bone 
geometry and stacking sequence [0,/90,],. Details on the tension-tension fatigue tests are available in [16]. The data 
acquisition stopped after Nr = 100 000 load cycles characterized by a load frequency of 5 Hz, sinusoidal shape, 
maximum force F = 31 kN and load ratio R ~ 0.14. The outer coupon dimensions are 152.4 mm X 254 mm [width 
x height], and the properties of the plies made of Torayca T700G are expressed in Table 1. 


Table 1. Ply properties. 





Young ’s modulus E, [GPa] 127.55 
Transverse elastic modulus E> [GPa] 8.41 
In-plane Poisson's ratio V12 [-] 0.309 
In-plane shear modulus Gy2_~—- [GPa] 6.2 
Out-of-plane shear modulus G23 ~~ [GPa] 2.82 
Thickness t [mm] 0.152 





A series of X-ray images collected during the test allowed the estimation of the inner damage (matrix crack density 
and delamination), and a strain gauge rosette on the outer surface was used to record the longitudinal strain, Figure 
2 





(a) (b) 
Figure 2. Notched cross-ply coupon under test (a), and X-ray image (b). The light gray region close to the notch is the 
delaminated region, while the horizontal, light gray lines are the matrix cracks that span the width of the 90° sub-laminate. 


A post-processing phase of the X-rays was executed to measure the amount of matrix crack density, delamination, 
and the reduction of the stiffness. The work in [3] summarizes the procedure to extract the amount of damage from 
the X-rays. Given the first damage assessment, the algorithm provides an estimation of the RUL which is 
systematically updated once a new measurement becomes available. The number of samples used to run the 
algorithm is N, = 15000. The initial values of the noise variances and the model parameter vector are reported in 
Table 2. The model parameters have been initialized using the data coming from other coupons. The critical 
threshold is set equal to the amount of damage observed at Nr = 100 000 load cycles, which is X¢x = [Pcr, Der, 
Scr]? = [422, 0.0229, 0.88]". 


Table 2. Initialization of random noises and model parameter vector. 
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The filtered estimate of the damage progression is shown in Figure 3. Every time a new observation becomes 
available, the updating of the weights wi? consents to estimate the posterior distribution of the system’s state. The 
expected values and the confidence bands of the filtered estimates are representative of the goodness of the 
algorithm in monitoring the multiple, concurrent damage mechanisms. The matrix crack density seems to be well- 
estimated by the algorithm; the expected value approaches the observed matrix crack density as time passes by. 
Delamination is slightly underestimated by the algorithm for most of the fatigue life, while the estimated trend of the 
normalized stiffness has very narrow confidence bands, insomuch as most of the experimental observations fall 
outside the 99% of the o-band. Though, the observations of the stiffness degradation appear noisy, and the estimated 
trend seems to fall in between the observations. Then, the algorithm seems able to filter out disturbances and 
concentrates the samples in between the noisy observations correctly. 
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Figure 3. Posterior estimation of the damage growth against load cycles; matrix crack density (a), delamination (b) and 
normalized stiffness (c). 


After every posterior estimation of the damage state, the samples are further projected in the future using the 
prediction equation (20) to predict the RUL. The sequential estimation of the RUL is reported in Figure 4. The RUL 
prediction is already close to the true RUL at the very beginning of the fatigue life, after few load cycles. Though, 
the wide confidence band of the prediction suggests that the information is characterized by large uncertainty. Then, 
the confidence band shrinks over time, but the average RUL moves away from the true value (between 10 000 and 
40 000 load cycles). This implies that the value of the model parameters is sensibly changing during this stage to 
better fit the available data. After N = 40 000 load cycles, the expected RUL converges to the correct value. The 
confidence band also keeps shrinking around the average RUL, thus improving the prediction performance. 
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Figure 4. RUL prediction. The two additional lines indicating the band RUL + 10% of the end of life help in evaluating the 
goodness of the prediction. 


5. Conclusions 


This work proposes a model-based Bayesian framework for composite laminates exhibiting concurrent damage 
mechanisms. An energy-based model able to estimate the SERR and the longitudinal elastic modulus constitutes the 
core of the fatigue damage accumulation model, and the latter is embedded into a Bayesian framework based on 
particle filtering. The approach has been already presented in the recent literature for monitoring the matrix crack 
density, but the framework proposed here enables the real-time monitoring and prediction of coexisting damage 
modes, which interact with one another and their combined effect influences the remaining life of the structure. The 


methodology has been successfully applied to CFRP damage progression data obtained through tension-tension 
fatigue experiments. 


The filtered estimation of the damage extent emphasizes that the algorithm slightly underestimates the observed 
delamination growth. The posterior expectation of the matrix crack density is in line with the observed p, and the 
estimated stiffness degradation appears well-centered on the noisy observations made with the strain gauge. High 
fluctuations of the RUL predictions characterize the early-mid stage of the fatigue life, but the expected value 
successfully converges to the true RUL after 50 000 load cycles and remains close to the target RUL until the end of 


the test. 


Future research should include the application of the prognostic method to other fatigue damage progression data. 
This is an importantstep to verify algorithm performance and generalize the validity of the approach. The number of 
samples N, has been chosen with a trial & error approach, and it is likely not the overall best choice to explore the 
state-space correctly. The use of refined methods to select N; can improve the algorithm performance or save 
computational time. In line with this idea, real-time, non-parametric methods to select (or update) the variance of the 
artificial dynamics perturbation can enhance the model parameter updating. Also, the use of refined sub-algorithms 
for the real-time parameter updating, like the kernel smoothing method may speed-up the convergence of the RUL 
prediction. Eventually, the definition of the end-of-life of CFRPs according to the asymptotic damage behavior 
should be further investigated. As a matter of fact, the horizontal asymptote representing the critical damage can 
sensibly change with the coupon because of the natural, inherent variability of the material. Therefore, a general and 
unified guideline to define the RUL of a composite laminate would be desirable. 
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